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Abstract 

Thresholding iterative methods are recently successfully applied to image deblurring prob¬ 
lems. In this paper, we investigate the modified linearized Bregman algorithm (MLBA) used 
in image deblurring problems, with a proper treatment of the boundary artifacts. We con¬ 
sider two standard approaches: the imposition of boundary conditions and the use of the 
rectangular blurring matrix. 

The fast convergence of the MLBA depends on a regularizing preconditioner that could be 
computationally expensive and hence it is usually chosen as a block circulant circulant block 
(BCCB) matrix, diagonalized by discrete Fourier transform. We show that the standard 
approach based on the BCCB preconditioner may provide low quality restored images and 
we propose different preconditioning strategies, that improve the quality of the restoration 
and save some computational cost at the same time. Motivated by a recent nonstationary 
preconditioned iteration, we propose a new algorithm that combines such method with the 
MLBA. We prove that it is a regularizing and convergent method. A variant with a stationary 
preconditioner is also considered. Finally, a large number of numerical experiments shows 
that our methods provide accurate and fast restorations, when compared with the state of 
the art. 


1 Introduction 

Image deblurring is the process of reconstructing an approximation of an image from blurred 
and noisy measurements. By assuming that the point spread function (PSF) is spatially- 
invariant, the observed image g{x,y) is related to the true image f{x,y) via the integral 
equation 

+ 00 +00 

9{x,y)= J j h{x-x',y-y')f{x\y') dx'dy'+ r]{x,y), (a;, y) G II C (1) 
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where r]{x, y) is the noise. 

By collocation of the previous integral equation on a uniform grid, we obtain the grayscale 
images of the observed image, of the true image, and of the PSF, denoted by G, F, and H, 
respectively. Since collected images are available only in a finite region, the field of view 
(FOV), the measured intensities near the boundary are affected by data outside the FOV. 
Given an n x n observed image G (for the sake of simplicity we assume square images), and 
a p X p PSF with p < n, then F is m x m with m = n + p — 1. Denoting by g and f 
the stack ordered vectors corresponding to G and F, the discretization of © leads to the 
under-determined linear system 

g = Ai + T], (2) 

where the matrix A is of size x When Imposing proper Boundary Conditions (BCs), 
the image A becomes square x and in some cases, depending on the BCs and the 
symmetry of the PSF, it can be diagonalized by discrete trigonometric transforms. For 
example, the matrix A is block circulant circulant block (BCCB) and it is diagonalizzable 
by Discrete Fourier Transform (DFT), when periodic BCs are imposed. 

Due to the ill-posedness of ([T]), A is severely ill-conditioned and may be singular. In such 
case, linear systems of equations ([H are commonly referred to as linear discrete ill-posed 
problems [53]. Therefore a good approximation of f cannot be obtained from the algebraic 
solution (e.g., the least-square solution) of ([5]), but regularization methods are required. The 
basic idea of regularization is to replace the original ill-conditioned problem with a nearby 
well-conditioned problem, whose solution approximates the true solution. One of the popular 
regularization techniques is the Tikhonov regularization and it amounts in solving 

mm{||Af-g||^-hAi||f||^}, (3) 

where || • ||p denotes the vector p-norm, p > 1, and p > 0 is a regularization parameter to 
be chosen. Hereafter, we use || • || = 11-112 to denote the t' 2 -uorm. The first term in 
is usually refereed as fidelity term and the second as regularization term. This approach 
is computationally attractive, since it leads to a linear problem and indeed several efficient 
methods have been developed for computing its solution and for estimating p [24] . On the 
other hand, the edges of restored image are usually over-smoothed. To overcome this un¬ 
pleasant property, nonlinear strategies have been employed, like total variation (TV) [35] 
and thresholding iterative methods [nils]. Anyway, several nonlinear regularization meth¬ 
ods have an inner step that apply a least-square regularization and hence can benefit from 
strategies previously developed for such simpler model, as we will show in the following. 

In this paper we consider a regularization strategy based on wavelet decomposition that 
has been recently largely investigated [SlEIIIlEsl HH Us]. This approach is motivated 
by the fact that most real images usually have sparse approximations under some wavelet 
basis. In particular, in this paper we consider the tight frame systems previously used in 
ISEIE]. Solving @ in a tight frame domain, the redundancy of system leads to robust signal 
representation in which partial loss of the data can be tolerated without adverse effects. In 
order to obtain the sparse approximation, we minimize the weighted £i-norm of the tight 
frame coefficients. Let be a wavelet or tight-frame synthesis operator {W’^W = /), the 
wavelets or tight-frame coefficients of the original image f are x such that 

f = IT^x. (4) 

In the following, we will investigate the synthesis approach, but our proposal can be applied 
also to the analysis and to the balanced approach described in E1I3Z]. Reformulating the 
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deblurring problem ([2]) in terms of frame coefficients 

min{/r||x||i + ||xf :Abb^x = g}, (5) 

X 

a regularized solution of this problem can be obtained by the Bregman splitting [32] ■ As for 
the iterative soft-thresholding [iiiia for the unconstrained version of ([S]) and the Landweber 
method for the least-square solution of the Bregman splitting converges very slowly for 
image deblurring problems. Hence a preconditioning strategy is usually employed, obtaining 
the Modified Linearized Bregman Algorithm (MLBA) (5). The preconditioner is usually 
chosen as a BCCB approximation of -I- al) Of > 0, see (5] jS] [57] . which is the 

simplest regularized version of the inverse of AA^ which, also when theoretically available, 
cannot be computed due to the severe ill-conditioning of A. 

In this paper, we show that the BCCB preconditioner used in the literature leads often to 
poor restorations when the matrix A has the rectangular A x rn? structure or is obtained by 
imposing accurate BCs, like antireflective BCs [35] ■ Note that in real applications we have 
to take into account the boundary effects to obtain high quality restorations, otherwise the 
restored image is severely affected by ringing effects [IZIISI]. This topic has been recently 
investigated in connection with nonlinear models based on wavelets or TV in [35] mo, but, 
at our knowledge, this is the first time that it is considered in connection with the MLBA. 
In this context we propose and discuss other preconditioning strategies for the MLBA with 
the synthesis approach. Our preconditioners are inspired by the experience with least-square 
regularization where the regularization preconditioning is studied since a long time, see the 
seminal paper [35]. In particular a nonstationary preconditioned iteration suggested by m 
leads to a new algorithm that is no longer a Bregman iteration. 

We investigate the following two strategies to define accurate and computationally cheap 
preconditioners: 

(1) an approximation of the blurring operator in a small Krylov subspace; 

(2) a symmetrization of the original PSF H; 

The choice (1) is quite natural and already considered for similar problems (see e.g. [2]), 
but we will show that a properly chosen Krylov subspace of small size (say spanned by at 
most five vectors), with a proper choice of the initial guess, is usually enough to obtain a 
good approximation. The choice (2) can be very useful in many applications where the PSF 
is obtained experimentally by measurements and is a perturbation of a symmetric kernel. 
In this case the approximated quadrantally symmetric (i.e., symmetric with respect to each 
quadrant) PSF leads to a matrix diagonalizable by Discrete Cosine Transform (DCT). 

Following the idea to combine preconditioned regularizing iterative methods for least- 
square ill-posed problems with the MLBA, we propose a new algorithm based on the recent 
proposal in [ISIIIH]. The nonstationary preconditioner is defined by a parameter computed 
by solving a nonlinear problem with a computational cost of 0{n^). We observe that this 
method can be applied only with square matrices and so only when A is obtained by im¬ 
posing BCs. The new algorithm is no longer a Bregman iteration and we cannot apply the 
convergence analysis developed in [5] for the MLBA. Therefore, here we prove its convergence 
and its regularization character. Furthermore, when a good value for the parameter in the 
preconditioner is available, we provide a variant of the algorithm with a stationary precon¬ 
ditioner which can improve the quality of the restorations even if the previous convergence 
analysis does not hold any longer. 

A large number of numerical experiments in Section [5| shows that our proposals not only 
outperform the standard MLBA with BCCB preconditioning, but are also good competitors 
for other recent methods dealing with boundary artifacts proposed in Hi. 
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Besides, we mention that the two deblurring models based on the rectangular matrix A 
and the imposition of BCs, we have tested also a third strategy based on the enlargement of 
the domain to reduce the ringing effects like in [Ml EH], but, according to the results in [2], 
the quality of the restored images were not better than those obtained with the other two 
models, while the CPU time was higher. Hence, we do not discuss further this strategy here. 

The paper is organized as follows. In Section [21 we describe the structure of the blurring 
matrix A explaining how fast trigonometric transforms can be used in the computations 
both for the rectangular matrix and the square matrices arising from the imposition of 
classical BCs. Section [3] reviews briefly the synthesis approach and the MLBA for solving 
the corresponding minimization problem. In Section |4] we propose possible regularization 
preconditioners, combining accurate restorations and a low computational cost. A new 
algorithm is proposed in Section |S] combining the MLBA with the method in [TB]. Section [5] 
contains a large number of numerical experiments, comparing our proposal with some state 
of the art algorithms, for the restoration of images with unknown boundaries. Concluding 
remarks are provided in Section |7l 


2 The structure of the blurring matrix 

We count mainly three strategies in order to obtain both accurate and fast restorations with 
reduced boundary artifacts. In this paper, we just consider two of them: the use of the 
original rectangular matrix and the imposition of BCs. As mentioned in the Introduction, 
we do not consider the third strategy introduced in [34], since from one side it is equivalent 
to the reflective (or Neumann) BCs in the case of quadrantally PSF, cf. [20], and from the 
other side, according to several tests that we have performed and the numerical results in 
[2], it does not provide a better restoration than the other two strategies, while it usually 
requires a larger CPU time. 

In this section we describe the structure of the matrix A and how fast computations 
with such matrix, like matrix-vector product or least-square solutions, can be implemented. 
Firstly we introduce the rectangular x m? matrix and then the square v? x matrix 
obtained when imposing proper BCs. 

2.1 The rectangular matrix 

The fact that this matrix is not square prevents the use of Fast Fourier Transform (FFT). To 
cope with this difficulty, one can construct an m? x rn? blurring matrix Aug that is BCCB, 
and hence, the FFT can be used. Let M £ M" be a masking matrix which, when applied 
to a vector in , selects only the entries in the FOV, i.e., their rows are a subset of the 
rows of an identity matrix of order rn?. The rectangular blurring matrix can be written as 

A = MAbig. (6) 

Hence the matrix vector Ax can be easily computed by two bi-dimensional FFTs of order 
TO^, followed by a selection of the pixels inside the FOV. Similarly A^x = A^j^M^x and 
thus two FFTs are applied to the zero-padded version of x of size rri^. This approach was 
used in [2] and it is numerically equivalent to that adopted in [10] . 

We observe that the matrix M in ([6]) leads to a matrix A independent of the BCs used 
in the definition of Abig. Thus, when the PSF is quadrantally symmetric, we suggest to use 
of the DCT instead of the DFT (see the discussion on reflective BCs in the next subsection). 
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Table 1: Pad of the original image F obtained by imposing the classical BCs considered in |27j . 
with Fc = fliplr(F), F^. = flipud(F), and F^c = flipud(fliplr(F)), where fliplr(-) and flipud(-) are 
the MATLAB functions that perform the left-right and up-down flip, respectively. 


The structure of the matrix A and its representation in ([5]) allow fast computations of the 
matrix vector product with A and but they prevent the use of fast transforms for solving 
linear systems with polynomials of AA^ as coefficient matrix even when A is full-rank. 

2.2 Boundary conditions 

The BC approach forces a functional dependency between the elements of f external to the 
FOV and those internal to this area. This has the effect of extending F outside of the FOV 
without adding any unknowns to the associated image deblurring problem. Therefore, the 
matrix A can be written as an v? x square matrix, whose structure can be exploited by 
fast algorithms. If the BC model is not a good approximation of the real world outside the 
FOV, the reconstructed image can be severely affected by some unwanted artifacts near the 
boundary, called ringing effects [27]. 

The use of different BCs can be motivated from information on the true image and/or 
from the availability of fast transforms to diagonalize the matrix A within 0(n^log(n)) 
arithmetic operations. Indeed, the matrix-vector product can be always computed by the 
2D FFT, after a proper padding of the image to convolve (the resulting image is the inner 
n X n part of the convolution), cf. [33], while the availability of fast transforms to diagonalize 
the matrix A depends on the BCs. Anyway, the shift-invariant property of the blur leads to 
a matrix A that can be well approximated by a BCCB matrix C, which is diagonalized by 
DFT, because usually in the applications 

A-C = R + N, (7) 

where i? is a matrix of small rank and A is a matrix of small norm. More precisely, for any 
£ > 0 there is a constant > 0 independent of n and depending only on e and on the PSF, 
such that the splitting © holds with 

rank(i?) < Cg • n, || A|| < £, (8) 

where rank(i?) denotes the rank of R (see [H]). Note that v? is the size of the matrix A. 

In the following we recall common BCs that will be used in the numerical tests. For a 
detailed description of zero, periodic, and reflective, refer to [22], while for antireflective BC 
see the review paper [21] and the original proposal in [36]. 

Zero (i.e., Dirichlet) BCs assume that the object is zero outside of the FOV. That is, one 
assumes that F has been extracted from a larger array padded by zeros (see Tabl^. 
This is a good choice when the true image is mostly zero outside the FOV, as is the 
case for many astronomical or medical images with a black background. Unfortunately, 
these BCs have a bad effect on reconstructions of images that are nonzero outside the 
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border, leading to reconstructed image with severe “ringing” near the boundary. The 
corresponding matrix A has a block-Toeplitz-Toeplitz-block (BTTB) structure which 
is not diagonalizable by fast trigonometric transforms. 

Periodic BCs assume that the observed image repeats in all directions. More specifically, 
one assumes that F has been extracted from a larger array of the form in Table [1] The 
corresponding matrix A is BCCB and so is always diagonalized by 2D DFT. Clearly, if 
the true image is not periodic outside the FOV the reconstructed image will be affected 
by severe ringing effects. 

Reflective (i.e., Neumann or symmetric) BCs assume that outside the FOV the image is a 
mirror image of F in]- That is, one assumes that F has been extracted from a larger 
array symmetrically padded like in Table [H The matrix A has a block structure that 
combines Toeplitz and Hankel structures, but it can be easily diagonalized by the DCT, 
when the PSF is quadrantally symmetric [81) . 

AntiReflective BCs have a more elaborate definition, but have a simple motivation: the 
antisymmetric pad yields an extension that preserves the continuity of the normal 
derivative [86) . They are given by 


F(l-i,j) = 2F(l,j)-F(i + l,J), 
F(i,l-j) = 2F(i,l)-F(i,j + l), 
F(n + i,j) = 2F{n,j) - F{n - i,j), 
F{i, n + j)= 2F{i, n) - F{i, n - j), 


I < i < p, 1 < j < n; 
1 < i < n, 1 < j < p; 


I < i < p, 1 < j < n; 
1 < z < n, 1 < j < p; 


for the edges, while for the corners the more computationally attractive choice is to 
antireflect first in one direction and then in the other m- This yields 

F{1 -t,l-j) = 4F(1,1) - 2F(1, j + 1) - 2F{i + 1,1) + + 1, j + 1), 

F{1 — i,n + j) = 4F(1, n) — 2F{1, n — j) — 2F{i -|-1, n) -|- F{i -|- 1, n — j), 

F{n + i,l — j) = 4F(n, 1) — 2F{n,j + 1) — 2F{n — z, 1) -|- F(n — i,j + 1), 

F{n + i,n + j) = 4F(n, n) — 2F{n, n — j) — 2F{n — i,n) + F{n — i,n — j), 

for 1 < z, J < p. 

The structure of the matrix A is quite involved, but it can be diagonalized by the antire- 
flective transform, when the PSF is quadrantally symmetric [T]. Since A is not normal 
the antireflective transform is not unitary, but it can be represented as a modihcation 
of the discrete sine transform, formed by adding a uniform sampling of constant and 
linear functions to the eigenvector basis preserving an “almost” unitary behaviour 
Due to the structure of A, the application of could generate artifacts at the boundary 
[19]; consequently it was proposed to replace by a reblurring matrix A obtained 
by imposing the same BCs to the PSF rotated by 180° m- Note that in the case of 
zero and periodic BCs A = A^. Furthermore, the MATLAB Toolbox RestoreTools 
[30) (that we will use in the numerical results) implements the reblurring approach 
to overload the matrix-vector product with A^, see [Ej. Therefore, for the sake of 
notational simplicity and uniformity, in the following the symbol has to be intended 
as the reblurring matrix A' in the case of antireflective BCs. 

^MATLAB functions for working with antireflective BC (antireflective transform, antisymmetric pad, ecc.) 
can be download at http://scienze-como.uninsubria.it/mdonatelli/Software/software.html 
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The reflective BCs will not be considered in the numerical results in Section [6] since they 
have the same computational properties of the antireflective BCs, e.g., for quadrantally PSF 
simply replace the antireflective transform with the DCT, and usually provide restorations 
slightly worser or at least comparable with the antireflective BCs, as we have numerically 
observed according to similar results with other regularization strategies [TOl [221 m] • On the 
other hand, we have recalled also the reflective BCs to motivate a preconditioner based on 
the DCT, instead of the DFT, when the PSF is quadrantally symmetric. 


3 MLBA for the synthesis approach 

For the synthesis approach [21122] the coefficient matrix is 

K = (9) 

where < s, A is the blurring matrix and is a tight-frame or wavelet synthesis operator. 
Note that using tight-frames W^W = I but ^ I m- The use of tight-frames instead 

of wavelets is motivated by the fact that the redundancy of tight-frame systems leads to 
robust signal representations in which partial loss of the data can be tolerated, without 
adverse effects, see e.g. [5]. 

Denote by x the frame coefhcients of the image f according to ([4]). Let the nonlinear 
operator be defined component-wise as 

[S^(x)], = (10) 


with Sfi the soft-thresholding function 

= sgn(a:i) max{|xi| - fi, 0} . 


Note that for image deblurring problems the singular values of A, and so those of K, decay 
exponentially to zero and we cannot assume that K is surjective. Therefore, the deblurring 
problem can be reformulated in terms of the frame coefficients x as 


mm L||x||i + ^||xf : arg mm ||i^x - gf 


( 11 ) 


which is equivalent to ([2]) if iF is surjective. 
The following linearized Bregman iteration 


J = z" -I- (g — ifx"), 

\ x"+i = AS^(z"+i), 


( 12 ) 


where z° = x° = 0, was introduced in [32] to solve problem (ED and later applied to image 
deblurring in [2]. A detailed convergence analysis of the linearized Bregman iteration ED 
was given in [7] when K is surjective, but we report here Theorem 3.1 in [2] that does not 
require such assumption. 

Theorem 1 ([2])- Let K G -n? < s and let 0 < X < ■ Lhen the sequence 

{x^+i} generated by ED converges to the unique solution of ED- 

As observed in |7] the convergence speed of (iT2]) depends of the condition number of K 
which, as observed before, is very large for image deblurring and hence the method results 
to be very slow. To accelerate its convergence in the case of KK^ /, in [2] the authors 
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modified iteration (IT^ by replacing with itif, where denotes the pseudo-inverse of 
K. If K is surjective = K"'"{K since < s. For image deblurring problems they 

suggested to replace with a symmetric positive definite matrix P such that 


{KK^y = {AA^y. 

Then the MLBA for frame-based image deblurring becomes [5] 

J = z”-I-— ifx”), 

\ x"+i = AS^(z-+i), 


(13) 


(14) 


where z° = x° = 0 . 

Remark 2. The MLBA (1141) is the linearized Bregman iteration (1121) for the linear system 

which is a preconditioned version of original linear system Kx = g by the preconditioner 
P^/^. In fact, by replacing K and gin © by P^^^K and P^/^g, respectively, we obtain da. 

The previous remark is the key observation used in [S] to prove that the MLBA algorithm 
converges to a minimizer of 


min {mIIxIIi -f : x = arg min |1 ATx - g|||, 1 . (15) 

xGM I Za J 

Theorem 3 ([5]). Assume P is a symmetric positive definite matrix and let 0 < X < 
. Then the sequence {x"'+^} generated by the MLBA (fTH) converges to the unique 
solution of m- 

The standard choice for P is 

P = {KK'^+ al)-^ = {AA^+ al)-^. (16) 


In such case 

||a:^pa:|| < i 

and hence A = 1 is a good choice according to Theorem O When AA^ + al is not easily 
invertible other choices as P can be explored, but usually the quantity IliF^PArH becomes 
hard to estimate. Therefore, assuming that P is a good approximation of (AA^ -|- al)~^, we 
set A = 1 in the algorithm. Anyway, the validity of the assumption ||A'^PAr|| < 1 can be 
guaranteed choosing a large enough. 

In conclusion, we consider the following version of the MLBA for the synthesis approach 

r z"+i =z"-klFA^P(g-AW^x"), 

\ x"+i = S^(z"+i), 

stopped by the discrepancy principle as in [^, i.e., at the first iteration n = h > 0 such that 

||r"||< 7 - 5 <||r"||, n = 0 , 1 ,..., h - 1 , (18) 

where 7 > 1, 5 = ||r 7 ||, and r” = g — AVF^x" is the residual at the n-th iteration. Here 
z° = x° = 0 and we assume that the noise level S is explicitly known. 



4 On the choice of the preconditioner P 

In this section we explore possible choices of P {AA^ + al)~^ which are computationally 
attractive. Let A be the rectangular, anti-reflective or BTTB matrix depending on the 
chosen treatment of the boundary of the image, and let C be the BCCB obtained from 
the same PSF. Since the matrix P in Theorem [3] serves as a preconditioner to accelerate 
the convergence, in this section we describe some preconditioning strategies, in order to 
combine fast computations with accurate restorations achievable when P Ri {AA^ -|- al)~^. 
The first proposal in Section 14.11 is the classical approach already used in the literature, 
cf. [3 EJ |37] . The second proposal in Section 14.21 is an approximation strategy considered 
for similar methods, c.f. [5], but this is the first time that it is explored with MLBA. The 
third proposal is inspired by a similar approach used with numerical methods that require 
symmetric matrices, see [laisT]. 

4.1 BCCB preconditioner 

Let C be the matrix obtained imposing periodic BCs. As described in Section [2T2l the matrix 
C is diagonalizable by DFT. Hence, the matrix-vector product with the matrix 

P = {CC^ + al)-^ 

can be efficiently computed by FFT and its use was previously proposed in [3. 

Algorithm 1. 

r z"+i =z"-tlFA^(C'C^-ba/)-i(g-41F^x’"), 

\ = S^(z"+i). 

We suggest to replace the DFT with the DCT, in the case of a quadrantally symmetric 
PSF like the Gaussian blur. The latter choice not only is motivated by computational 
considerations, since the complex DFT is replaced by a real transform (DCT), but also by 
the quality of the computed approximation. Indeed, using the DCT the matrix C can be seen 
as an approximation of A by imposing reflective BCs instead of periodic BCs, which results 
usually in a better approximation and so provides better restorations. The same expedient 
will be used for the following preconditioners as well. 

Note that ||(C'C'’^-|-a/)“^|| < ||(AA’^-|-a/)“^|| is a sufficient condition to apply the The¬ 
orem [3] with A = 1. This assumption could be hard to be satisfied in practice. Nevertheless, 
it is expected that 

\\K'^{CC'^ + aI)-^K\\<l (20) 

if a is large enough. 

In the literature regarding preconditioning of Toeplitz matrices by circulant matrices, 
several strategies have been proposed to compute the matrix C, cf. [9]. In this paper 
we simply consider the matrix obtained imposing periodic BCs, which corresponds to the 
natural Strang preconditioner, since we have not observed numerical differences, when using 
other strategies like the optimal Frobenius norm approximation preconditioner Roughly 
speaking, this follows from the fact that the entries of A depend on the value of the pixels 
of the PSF according to the shift invariance structure. In particular the central coefficient 
of the PSF belongs to the main diagonal of A and the pixels near the center of the PSF 
belongs to the central diagonals of the central blocks. Finally, the PSF is almost centered 
in the middle of a n x n image and hence every pixel is distant at most n/2 pixels in every 
direction (usually much less due to the compact support). Hence the block band and the 
band of each block of A are at most n/2. 
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4.2 Krylov subspace approximation 

The preconditioner in the previous section is essentially defined as an approximation of the 
operator A in the Fourier domain. Another strategy, useful also for more general matrices, is 
to employ orthogonal or oblique projections into subspaces of small dimension. A common 
choice is a proper Krylov subspace. 

The matrix vector product t" = (AA^ + can be computed solving the linear 

system 

(AA'^ + a/)t" =r", (21) 

whose solution can be approximated by few iterations of conjugate gradient (CG) since 
AA^ + al is symmetric and positive definite. One or few steps of CG to approximate the 
vector (AA^ + a/)“^r" is a common strategy, see e.g. [5]. Here we explore the use of a good 
preconditioner associated with a proper choice of the initial guess and the stopping criteria. 
We solve m by preconditioned CG (PCG) with preconditioner the matrix (CC^ + al) ^ 
introduced in Section 14.11 This is equivalent to solve the linear system 

{CC^ + aI)-^/‘^{AA^ + aI){CC^ + = {CC^ + (22) 

with y" = {CC^ + a/)i/2t". 

The Krylov subspace of size j generated by the matrix B and the vector v is defined by 
/Cj(il, v) = spanjv, Hv,..., jG N. 


We denote by the vector that minimizes the energy norm of the error of the linear system 
(|^ into the Krylov subspace 

■= ((CC^ + aI)-^/‘^{AA^ + aI)iCC^ + {CC^ + . 


Therefore, defining 

the following algorithm can be sketched. 


Algorithm 2. 


r = z" + WA^t^ , 

\ x"+i = S^(z"+1). 


(23) 


Of course a large Pn is not practical and also the convergence of the Algorithm 2 could fail 
if is not a good approximation of t” in ([71]) . However, in practice Pn can be taken very 
small and the PCG converges very rapidly assuring also the convergence of the Algorithm 2 
as numerically confirmed by the results in Section [51 This follows from discussion at the 
beginning of Section 12.21 and the well-conditioning of the coefficient matrix of the linear 
system (ED). Indeed, all the eigenvalues of AA^-|-a/ are in [a, c], with c constant independent 
of n and usually c = 1 -|- a, because A arises from the discretization of ([T|) and, thanks to the 
physical properties of the PSF (nonnegative entries and sum of all pixels equal to one), its 
largest singular value is bounded by one. If follows that the BCCB preconditioner CC"'" + al 
is very effective since the spectrum of {CC'^+aI)~^^'^{AA^ + aI){CC'^ + is clustered 

at 1, with 0{n) outliers according to equations d?]) and (|8]), while is the total number of 
eigenvalues (see [25l [T0] b 

Moreover, we observe that if the Algorithm 2 is converging then, in the noise free case, 
the residual is going to zero (otherwise stagnates around <5) and hence also the solution of 
the linear system (EH) approaches the zero vector. This has two interesting consequences. 
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First, a good initial guess for the PCG is the zero vector since it is a good approximation 
of the solution of the linear system at least for n large enough. Second, the size of the 
Krylov subspace should decreases when n increases reaching the same fixed accuracy in 
the approximation of t" (see the following discussion on /3„). 

Note that the computation of requires /3n matrix-vector products with AA^ + al and 
with CC'^ + al. Nevertheless, according to the previous discussion, a small /3„, e.g., /3„ < 5, 
is enough. In the numerical results in Section |6] we fix 


/3„ = min{5,/3toi}, (24) 

where /3toi is the number of PCG iterations required for reaching the tolerance 10“^ in 
terms of the norm of the relative residual in the linear system (1221) . We observe that in our 
numerical results, /3„ decreases quickly obtaining /3„ = 1 for all n > h, with h small. 

Finally, we note that the preconditioner P obtained by the PCG approximation is not 
stationary and changes at each iteration of the MLBA. Therefore, Theorem [3] cannot be 
applied. Nevertheless, Algorithm 2 with the condition (IMl) has been convergent in all our 
numerical experiments confirming that is a good approximation of t", at least for n large 
enough. 

4.3 Preconditioning by symmetrization 

The preconditioner in Section 14.11 is related to a different boundary model, namely periodic 
BCs, but the deblurring problem and in particular the PSF are the same. Unfortunately, 
periodic BCs lead to poor restorations for generic images and for more accurate models, like 
reflective or antireflective BCs, the matrix A cannot be diagonalized by fast trigonometric 
transforms when the PSF is not quadrantally symmetric. Therefore, in this section we 
use a different strategy: the preconditioner is defined by a different PSF that leads to fast 
computations with an accurate deblurring model. 

We consider a simple implementation of this strategy that can be useful when the PSF 
is experimentally measured. Indeed, in some applications, the PSF is nonsymmetric even if 
it is just a numerical perturbation of a Gaussian-like blur, cf. Example 1 and |25j . Recall¬ 
ing that for the reflective and antireflective BCs fast transforms (cosine and antireflective, 
respectivelly) can be implemented only in the quadrantally symmetric case, a quadrantally 
symmetric PSF H can be obtained from the original PSF H by dehning 

uf H{i,j) + H{-i,j) + H{i,-j) + . . 

= - - -, i,j = l,...,n. 

Note that H is the optimal Frobenius norm approximation of H in the set of quadrantally 
symmetric PSFs, see [81]. Therefore, we consider the matrix Q obtained imposing reflective 
BC to H when A is the BTTB or the rectangular matrix, while for A antireflective, Q is 
defined imposing antireflective BCs as well. In this way 

P= (QQ^ + a/)-i 

can be diagonalized by DCT or by antireflective transform and the MLBA becomes 

Algorithm 3. 


r z”+i = z" -h WA^{QQ'^ + al)-\g - AIU'^x"), 
\ X-+1 = S^(z-+i). 
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In analogy to Algorithm 1, \\{QQ'^ + al) ^|| < ||(AA^ + al) ^|| is a sufficient condition 
to apply Theorem [3] with A = 1. It is expected that 

\\K^{QQ^ + aI)-^K\\ < 1 


if a is large enough. 


5 Approximated Tikhonov regularization instead of pre¬ 
conditioning 

In this section we propose an approach to approximate different from the use of the 
matrix P in p3p as suggested in . Motivated by a very recent preconditioning proposal in 
[13 [IH], we replace the whole matrix by a regularized approximation obtained by C. 

In [15] the authors suggest to solve the preconditioned linear system 

ZAi = Zg, 

where Z is a regularized approximation of K\ by a Van Cittert iteration |39j . instead 
to solve a preconditioned Landweber iteration. Unfortunately, ZA is not symmetric and 
the convergence analysis, based on the complex eigenvalues of ZA, is hard to generalize. 
Differently, the nonstationary preconditioned iteration proposed in |18j results in a similar 
iteration, but an elegant convergence analysis is provided under a minor approximation 
assumption. 

For P = {AAJ" + al)~^, the correction term K^P{g — ATx") in the MLBA (fTdll can be 
seen as the Tikhonov solution, with parameter a, of the error equation. In detail, in 

the iteration o can be rewritten as 


^n+l 


(26) 


where 


p" = K'^P{g - ATx”) 

= K^{KK^ + al)-\g - Kyp) 

= (K'^K + aI)-^K'^{g - Kyp), 

since (AA^ + al)-^ = {KK^ + al)-^ and K^{KK^ + al)-^ = [K'^K + ar)-^K'^. Note 
that the correction p" is the solution of the Tikhonov problem 

mmlllATp-r^f+ a||pf}, 

which is a regularized approximation of the error equation 

iFe" = r" (27) 

in the noise free case, i.e., 5 = 0, where e” = x —x” denotes the error at the current iteration. 

In real applications S 0 and so equation (1771) is (only) correct up to the perturbation in 
the data. Taking this into account, one may as well consider instead of the error equation (1271) 
the “model equation” 

Le" = r" , (28) 


12 


where L is an approximation of K, possibly tolerating a slightly larger misfit. Solving (l28l) 
by means of Tikhonov regularization, we find 

p" = (L^L + a/)-iL^r" 

= WC'^iCC'^ + a/)-ir", 

where we have choosen 

L = CW^. (29) 

Using p" in (I26p to replace p” we obtain a new algorithm. 

Algorithm 4. 

r z"+i =z" + lUC''^(CC''^ + a/)-Hg- 
\ x"+i = S^(z"+1). 

As before, the matrix C is chosen as a BCCB in general or diagonalizable by DOT (i.e., 
the reflective BC matrix), when the PSF is quadrantally symmetric. Unfortunately, this 
preconditioning strategy cannot be applied to the rectangular matrix approach because, in 
such case, C should have the same size of A, but this condition prevents the possibility of 
computing p" by fast trigonometric transforms. 

Remark 4. The iteration (1301) uses the preconditioned linear system 

WC^{CC^ + a/)-i Ax = 1UC7'^(CC''^ + a/)“^g, (31) 

to update an aproximation inspired by the linearized Bregman iteration but without 

resorting to the normal equations. 

Clearly, Algorithm 4 is no longer a MLBA and so a different convergence analysis is re¬ 
quired. Unfortunately, classical results for convex optimization cannot be applied since the 
coefficient matrix 1UC^(C'C'^ + aI)~^K in (I5T|) is not symmetric positive dehnite. An alter¬ 
native convergence proof could be very hard because also the complex analysis convergence 
in [H] cannot be easily combined with the Bregman splitting and soft-thresholding. 

Therefore, accordingly to m. we consider a nonstationary choice of a that allows to 
provide a convergence analysis of the resulting algorithm and avoid the a-priori choice of a. 
On the other hand, if a good estimation of a is available, then Algorithm 4 can provides 
better restorations and hence it is also considered in the numerical results in Section [HI 

Assumption 1. Let A,C G and W G < s, such that 

||(C-A)v|| <p||Av||, VveR"', (32a) 

and 

||ClU^(u-5^(u))|| <p(5, VueR®, (32b) 

with a fixed 0 < p < 1/2, where S = ||? 7 || is the noise level. 

The assumption (I32al) is the same spectral equivalence required in [18] . Let L be defined 
in dMl), then equation (|32al) translates into 

|1(L-A)u|| <p||Au||, VueR®. (33) 

Instead, the assumption (I32bl) was not present in |18j and it is equivalent to consider the 
soft-threshold parameter p as a continuous function with respect to the noise level 6, i.e., 
pL = p{6), and such that p{6) —^ 0 as d —>■ 0. This is a common request in many soft- 
thresholding based methods, see for instance Theorem 4.1 in [T3|. Nevertheless, in this paper 
we will not concentrate on p and will not give any specific i5-dependent rule to compute it. 


13 







Algorithm. 4—NS. Let tP he given and set r° = g — A5'^(z°). Choose t = P 

from (I32a|) . and fix q G ( 2 / 9 ,1). 

While ||r”|| > tS, let Tn = ||r"||/iJ and let a„ be such that 

an\\{CC'^ + a„/)“^r"|| = g„||r”||, Qn = max{g, 2 p+ (1 + p)/t„}, 

compute 

I z"+i = z" + WC^iCC^ + a„I)-\g - 
\ x"+i = S^(z"+i). 

Note that the iteration (I34bl) is the same of Algorithm 4 where a nonstationary a is 
chosen at every iteration according to (I34a|) . In Corollary [ 8 ] we will prove that, if (5 > 0, then 
Algorithm 4-NS will terminate after n = ns > 0 iterations with 

||r”'’|| < rh < ||r"'||, n = 0,1, • • • , - 1, (35) 

which is the discrepancy principle (1181) with n = ns and 7 = r = (1 + 2 /o)/(l — 2p). 

The parameter q in Algorithm 4, like in |18| . is meant as a safeguard to prevent that the 
residual decreases too rapidly. Our theoretical results do not utilize this parameter. 

Remark 5. It is not difficult to see that there is a unique positive parameter a„ that satisfies 
p4al) . This parameter can be computed with a few step of an appropriate Newton scheme 
Accordingly, parameter an, and therefore Algorithm are well defined. 

We define 

h" = L^(LL^ + a„/)-i(g - KSfiz^)), (36) 

such that (I34bl) can be compactly rewritten as 

z"+i = z" + h", 
x”+i = S^(z’^+i). 

For the purpose of the subsequent convergence and regularization results, when <5 > 0, 
even if it will be always the case, we will highlight by the subscript S (for instance {x^}) 
the sequences generated by Algorithm 4-NS, starting from initial data = Ai + t] affected 
by noise, whereas we avoid the subscript (for instance {x"}) for the sequences generated 
starting from exact initial data g = Af, i.e., S — 0 . 

For the following analysis, instead of working with the error = x — x^, it is useful to 
consider the partial error with respect to z^, namely 

e? = x-z?. (37) 

Proposition 6 . Assume that the assumptions (1321) are satisfied for some 0 < p < 1/2. If 
Ikill > and we define Tn = ||r^||/5, then it follows that 

||r^ - Le^ll < ||r?|| < (1 - p)||r?||, (38) 

where e" is defined in 

Proof. In the free noise case we have g = Ax. As a consequence 

r? - Le? = g^ - Ax^ - L(x - z?) + Lx? - LSfiz^) 

= g'-g+(A-L)e? + L(z?-5^(z?)). 


(34a) 

(34b) 
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Using now assumptions (l32l) . in particular (l33)) . and ||g'^ — g|| < <5, we derive the following 
estimate 

||r? - Le^ll < ||g^ - g|| + \\{K- L)e^\\ + ||L(z? - ^^(z^))|| 

< l|g^-g|| +Pll^e^ll +PS 

<||g^-g||+p(||r?|| + ||g^-g||+d) 

< (1 + 2p)S + pllr^||. 

The first inequality in (IMl) now follows from the hypothesis 6 = ||r^||/r„. The second 
inequality follows from p + < p + ■ □ 

We are going to show that the sequence {x^} approaches x as d —>■ 0. The proof combines 
Proposition |6] with suitable modifications of the results in [18]. 

Proposition 7. Let ef he defined in (l37)) . If the assumptions (l32l) are satisfied, then ||e^|| 
of Algorithm f-NS decreases monotonieally for n = 0,1,... ,ns — 1. In particular, we deduee 

||e?f - + a„)-^r^||||r?||. (39) 

Proof. We have 

lle^f-ller^f = 2(e?,h")-||h'^f 

= 2{Le2, (CC^ + aniy^r^) - {r^,CC^{CC^ + a„/)-"r^) 

= 2(r^, {CC^ + a„/)-ir?) - {r'f,CC^{CC^ + a^ir^v^s) 

- 2(r? - Le^, {CC^ + a„J)-ir^) 

> 2(r^, + a„/)-ir?) - 2(r^, CC^{CC^ + a^ir^v^f) 

-2(r?-Le?,(CC^+a„/)-ir?) 

= 2a„(r^, {CC^ + a^ir^v^f) - 2{v'f - Le^, [CC^ + a„/)-ir?) 

> 2a„(r^, {CC^ + a„/)-2r^) - 2||r^ - Le^||||(CC^ + an/)-'r^|| 

= 2||(CC^ + a„/)-ir?|| (||a„(CC^ + an/)-'r^|| - ||r? - Le^) 

>2||(CC^ + a„/)-ir?|| (^g„||r^|| - ||r?||^ 

>^^^\\{CC^ + ar.I)-^rmrn, 

where the relevant inequalities are a consequence of equation (I34ap and Proposition [51 The 
last inequality follows from (I34ap and > r = (1 + 2p)/(l — 2p) for ||r^|| > t6. □ 


Corollary 8. Under the assumptions (I32p . there holds 

R 2 ' ^■5- 1 raa - 1 

l|e^ll> E ll(CC^ + ««/)-'r?||||r^||>c^ ||r^f 

^ n—0 n—0 

for some constant c > 0, depending only on p and q in (I34al) . 


(40) 
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Proof. The following proof is almost the same as Corollary 3 in [18], but we include it to 
make the paper self contained. 

The first inequality follows by taking the sum of the quantities in ([3^ from n = 0 up to 
n = ns — 1. 

For the second inequality, note that for every a > and every a G cr(C') C [0, |lC|p], 

with (j{C) being the spectrum of C, we have 


a a 

- ||C||2 + a 


(l + ||C||Va)-i >g„. 


and hence, 

a\\iCC^ + aI)-h^\\>q4rn, 

as ||r^|| > 0 for n < na- This implies that a„ in (I34al) satisfies 0 < «„ < thus 


||(CC^ + a„J)-ir^|| = ^||r?||> 


(1 Qn) I 

||CP ' 


Now, according to the choice of parameters in Algorithm 4-NS, we deduce 
1- Qn = min{l - q, 1 - 2p - (1 + p)/t„}, and 


1 - 2p- (1 + p)/t„ 


1 + 1 p ^ 1 + ‘Ip \ p p 

r Tn r T T 


Therefore, there exists c > 0, depending only on p and q such that l — qn> cUCp 
and 


||(CC^ + a„J)-ir?||>c 


8p2 


1 + 2p 


-1 


for n = 0,1, • • • ,715 — 1. 


-1 


Now the second inequality follows immediately. 


□ 


From the outer inequality of (l40|) it can be seen that the sum of the squares of the residual 
norms is bounded, and hence, if d > 0, there must be a first integer ns < oo such that (13511 
is fulfilled, i.e., Algorithm 4-NS terminates after finitely many iterations. 

Remark 9. Recalling that the soft-threshold parameter p is taken as a continuous function 
with respect to the noise level 5 such that p(S) 0 as S —>■ 0, then the operator g i—>■ z" is 
continuous for every fixed n. 

In the next theorem we are going to give a convergence and regularity result. 

Theorem 10. Assume that z° is not a solution of the linear system 


g = AVF^x, 


(41) 


and that 5m is a sequence of positive real numbers such that 5m ^ 0 as m ^ oo. Then, 
the sequence generated by the discrepancy principle rule (1351) . converges as 

m ^ oo to the solution of m which is closest to in Euclidean norm. 

Proof. We are going to show convergence for the sequence and then the thesis 

will follow easily from the continuity of >S'^( 5 ) and Remark 19] i.e.. 


1* ?^(<5m) 

lim Xr 
m—>oo 


lim 
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The proof of the convergence for the sequence can be divided into two steps: at step 

one, we show the convergence in the free noise case <5 = 0. In particular, the sequence {z”} 
converges to a solution of (l4T|l that is the closest to z°. At the second step, we show that 
given a sequence of positive real numbers dm 0 as m —> c», then we get a corresponding 
sequence {z^ | converging as m —>• cx). 

Concerning the first step of the proof, we will not give details since it can be just copied 
from [18] [Theorem 4]. Indeed, if d = 0, from Remark |9| it follows that = r", and the 
sequence {z"} coincides with the one generated by algorithm 1 in [THj. We just say that the 
main ingredients are the convergence of the sequence ||e"|| granted by Proposition Hand the 
convergence to 0 of the sequence ||r"||||(C'C'^ + Q!„/)“^r”||, since general term of a converging 
series from CorollaryjSj Moreover, in the free noise case the sequence {z"} will not stop, i.e., 
n —> c», since the discrepancy principle will not be satisfied by any n, in particular ns ^ oo 
for d —>■ 0. 

Hence, let x be the converging point of the sequence {z"} and let dm > 0 be a sequence 
of positive real numbers converging to 0. For every dm, let n = n(Sm) be the first positive 
integer such that (l35)) is satisfied, whose existence is granted by Corollary |8l and let 
be the corresponding sequence. For every fixed e > 0, there exists n = n(e) such that 

||x — z"|| < e/2 for every n > n(e), (42) 

and there exists d = d(e) for which 

||z" — z^ll < e/2 for every 0 < d < d, (43) 

due to the continuity of the operator g z" for every fixed n, see Remark H Therefore, 
let us choose m = m(e) large enough such that 5m < 5 and such that n(5m) > n for every 
m >rn. Such m does exists since dm —5" 0 and —>■ c» for d —0. Hence, for every m > Wi, 
we have 

<iidi 

= l|x-4„ll 

<||x-zl + ||z"-z"^||<e, 

where the first inequality comes from Proposition |7| and the last one from (l42l) and (l43l) . □ 

6 Numerical results 

In this section, we will show the numerical results for image deblurring using our proposed 
Algorithms 1-4. We compare them with some available deblurring algorithms, which im¬ 
plement a proper treatment of the boundary artifacts. In particular we consider two of the 
algorithms proposed in |5|, namely FA-MD for the Frame-based analysis model and TV-MD 
for the Total Variation model, and the Algorithm |3] called here FTVd since, in the case of 
nonsymmetric PSF, it reduces to an implementation of the algorithm in m with the trick 
described in . The codes of the previous algorithms are available at the web-page of the 
authors and we use the default parameters and stop conditions. The regularization parameter 
is chosen by hand in order to provide the best restoration (see the following discussion). 

Our tests were done by using MATLAB 7.11.0 (R2010b) with floating-point precision 
about 2.22 • 10“^® on a Lenovo laptop with Intel(R) Core(TM) i2 CPU 2.20 GHz and 2 
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GB memory. Assuming that the noise level is available or easily estimated, we stop all 
Algorithms 1-4 using the discrepancy principle ([T^ with 7 = 10“^^. Algorithm 4-NS is 
stopped according to the modified discrepancy principle (IMI) . Moreover, for the Algorithm 4- 
NS we set 

q = 0.5 and p = 10“^. 

Therefore, q and p do not need to be estimated. 

The accuracy of the solution is measured by the PSNR value, which is defined as 


PSNR = 201ogio 


255 • n 

iif~fir 


with f and f being the original and the restored images in the FOV, respectively. The initial 
guess of each algorithm is set to be the zero vector. 

To estimate p and a, since they are mutually dependent and they are related to the 
preconditioner, we fix a possible p (usually the results are not very sensible varying p ii a 
is properly chosen) and then the optimum a, which gives the largest PSNR, is chosen by 
trial and error. Possible strategies to estimate a will be investigated in future works. Only 
for Algorithm4-NS we pay a slightly more attention in the choice of p since this is the only 
parameter of the method. Similarly, for all the other methods considered for comparison, 
the regularization parameter is chosen by trial and error, as the one leading to the largest 
PSNR. 

We take only the more appropriate BCs for each example. In particular, if the image 
has a black background, like in astronomical imaging, we consider zero BCs, while when 
the image is a generic picture we use antirefiective BCs. In the following the “Algorithm 
x" is denoted by “Alg-BCcc” and “Alg-Rectx”, when A is obtained imposing BCs or is the 
rectangular matrix, respectively. We recall that Algorithm 4 is available only for the BC 
approach. 


6.1 Linear B-spline framelets 

The tight-frame used in our tests is the piecewise linear B-spline framelets given in [5]. 
Namely, given the masks 


= i [1,2,1], 6i = ^[l, 0, -1], 62 = i[-l, 2, -1], 

we define the ID filters of size n x n hy imposing reflective BCs 


Bn = - 
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B2 = - 
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(a) true image 


(b) PSF image 


(c) observed image (d) Alg-BCl 


(e) Alg-BC3 



(f) Alg-Rect2 (g) Alg-Rect3 (h) Alg-BC4-NS (i) TV-MD (j) FTVd 

Figure 1: Example 1: true image, PSF, observed image, and restored images. 


The nine 2D filters are obtained by 


Bij = Bi 


B, 


hj = 0,1,2, 


where 0 denotes the tensor product operator. Finally, the corresponding tight-frame analysis 
operator is 


W = 


B< 


0.0 


So, 


B- 


2,2 


Throughout the experiments, the level of the framelet decomposition is 4 like in and 
the level of wavelet decomposition is the one used in FA-MD. 


6.2 Example 1: Saturn image 

The first example is 256 x 256 Saturn image in Figure [T] (a) while the astronomical PSF is 
taken from the “satellite” test problem in [33] Figure [T](b). We add a 1% of Gaussian white 
noise to obtain the observed image in Figure [T](c). We assume zero BCs. 

Note that Alg-BC3 and Alg-Rect3 use the DCT for the preconditioner, while Alg-BC4 
like Alg-BCl and Alg-Rectl use FFT, since the PSF is not quadrantally symmetric. 

Table [5] reports the PSNR and the CPU time for the different algorithms. Note that 
Algorithm 1 provides a poor and time consuming restoration. Moreover, it requires a larger 
value of the parameter a with respect to the algorithms 2-4, which is necessary to satisfy 
condition (1201) and assure the convergence. The algortihms 2 and 3 have the largest PSNR 
with reasonable CPU time, in particular Alg-Rect3 seems to be a good choice and Alg-BC3 
gives a comparable restoration in about half time. Algorithm 4 gives a slightly lower PSNR 
even if the computed restorations are better than Algorithm 1 and the other algorithms from 
the literature, keeping also a low CPU time. 















Algorithm 

PSNR 

Iter. 

CPU time(s) 

Regular, parameter 

Alg-BCl 

30.97 

322 

200.99 

a = 0.045 

Alg-BC2 

31.60 

9 

18.18 

a = 0.0004 

Alg-BC3 

31.56 

10 

7.07 

a = 0.0005 

Alg-Recti 

30.95 

493 

948.94 

a = 0.07 

Alg-Rect2 

31.62 

8 

22.20 

a = 0.0003 

Alg-Rect3 

31.61 

7 

13.50 

a = 0.0003 

Alg-BC4 

31.49 

29 

16.56 

a = 0.0018 

Alg-BC4-NS 

31.25 

15 

10.32 

/i = 6 

FA-MD 

30.87 


90.85 

A = 0.001 

TV-MD 

31.17 


47.61 

A = 0.01 

FTVd 

30.50 


1.75 

1/a = 0.0013 


Table 2: Example 1: PSNR, number of iterations, and CPU time in seconds for the best 
regularization parameter (maximum PSNR) reported in the last column. For our algorithms 
/r = 10 except for Alg-BC4-NS. 



(a) Alg-BCl (b) Alg-Rect3 (c) Alg-BC4 (d) FA-MD (e) TV-MD 


Figure 2: Example 1: residual image g — Af, where f is computed by different algorithms. 


The algorithms in [2] (FA-MD and TV-MD) in this example lead to a larger CPU time, 
while the FTVd is very fast but the computed restoration is the worst. Figure [T] shows the 
corresponding restored images. To test the quality of the restorations, Figure [5] shows the 
residual images defined as g — Af, where f is the restored image. 


6.3 Example 2: Galcixy image 

We consider another astronomical example with the 256 x 256 image in Figure [3] (a) cor¬ 
rupted by oblique Gaussian blur taken from the “GaussianBlur422” test problem in [S^, see 
Figure |3] (b). We add a 2% of Gaussian white noise to obtain the observed image in Fig¬ 
ure [3] (c). We impose zero BCs and the computational properties of the different algorithms 
are the same as in Example 1. 

Table [3] shows that Alg-BC4 is the best algorithm, since it obtains about the same PSNR 
of Algorithm 2, but with about 1/4 of the CPU time. The variant Alg-BC4-NS with a 
nonstationary choice of the preconditioner results to be very effective and comparable with 
the other algorithms based on the BC model avoiding the choice of parameter a. Differently, 
the rectangular approach gives a slightly lower PSNR with a larger CPU time. Concerning 
the other methods, the same observations reported for Example 1 still apply. Figure [3] shows 
some of the corresponding restored images. 
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(a) true image (b) PSF image (c) observed image (d) Alg-BCl (e) Alg-BC3 



(f) Alg-Rect3 (g) Alg-BC4 (h) Alg-BC4-NS (i) FA-MD (j) FTVd 

Figure 3: Example 2: true image, PSF, observed image, and restored images. 


Algorithm 

PSNR 

Iter. 

CPU time(s) 

Regular, parameter 

Alg-BCl 

25.02 

21 

16.57 

a = 0.04 

Alg-BC2 

25.07 

12 

22.88 

a = 0.008 

Alg-BC3 

25.05 

27 

17.58 

a = 0.02 

Alg-Rectl 

24.91 

25 

42.14 

a = 0.06 

Alg-Rect2 

24.98 

12 

28.21 

a = 0.007 

Alg-Rect3 

24.96 

21 

36.18 

a = 0.01 

Alg-BC4 

25.06 

12 

6.21 

a = 0.008 

Alg-BC4-NS 

25.01 

29 

17.10 

M = 4 

FA-MD 

24.50 


77.05 

A = 0.02 

TV-MD 

24.55 


68.91 

A = 0.09 

FTVd 

24.62 


1.51 

1/a = 0.027 


Table 3: Example 2; PSNR and CPU time for the best regularization parameter (maximum 
PSNR). For our algorithms // = 10 except for Alg-BC4~NS. 







(f) Alg-Rect3 



(g) Alg-BC4 


(h) Alg-BC4-NS 


(i) FA-MD 


(j) TV-MD 


Figure 4: Example 3: true image, PSF, observed image, and restored images. 


6.4 Example 3: Boat image 

To set up a scenario of unknown boundaries, the observed image of size 196 x 196 is obtained 
convolving the full (256 x 256) image using arbitrary BCs (periodic, for computational con¬ 
venience) and then keeping only the pixels in the FOV 196 x 196 (i.e., those not depending 
on the BCs). The FOV is denoted by a black box in the true image in Figure |4] (a). 1% of 
white Gaussian noise is added to the 196 x 196 blurred image. 

We impose antireflective BCs owing to the generic structure of that picture. Hence the 
preconditioner in Alg-BC3 is diagonalized by the antireflective transform, according to the 
structure of the matrix A. Alg-Rect3 uses the DCT as usual, while Alg-BC4 like Alg-BCl 
and Alg-Rectl use FFT since the PSF is not quadrantally symmetric. 

Table m shows the PSNR and the CPU time for the best restorations shown in Figure 21 
Note that our algorithms with the rectangular matrix are less effective than the antireflective 
BC approach, leading to a lower PSNR. To obtain reasonable restorations in the rectangular 
case we need a large and so we take a different /r for the two deblurring models (BC and 
rectangular matrix). The best algorithm results to be Alg-BC4, since it combines a good 
restoration with a low CPU time. 

6.5 Example 4: Cameraman with Gaussian blur 

In this example we consider the classical Cameraman image 256 x 256 distorted by a 31 x 31 
Gaussian blur with standard deviation 2.5 and a 2% of white Gaussian noise (see Figure [5|). 
The size of the observed image is 226 x 226 according to the support of the PSF. 

The PSF is quadrantally symmetric and hence, following our discussion in Section 14.11 
Algorithm 1 is implemented using the DGT instead of FFT. Clearly H = H and so Alg-Rect3 
reduces to Alg-Rectl (they are really the same algorithm). We impose antireflective BCs and 
hence Alg-BC3 is the standard MLB A (fT7)l with P = (AA^ + al)~^, but recalling that we 
are using the reblurring approach where is replaced by A. On the other hand Alg-BC2 
is no longer useful since the matrix-vector product with the matrix P = {AA^ + al)~^ can 
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PSNR Iter. CPU time(s) 


Regular, parameter 


^ = 200, a = 0.03 


30.17 13 3.67 

29.77 60 19.57 


a = 


Algorithm 
A 

a: 

_A_ 

Alg-Recti 

Alg-Rect2 

Alg-Rect3 

Alg-BC4 

Alg-BC4-NS 


g-BCl 

g-BC2 

e-BC3 


29.43 

97 

34.26 

30.11 

11 

17.21 

30.09 

10 

4.03 


a = 
a = 


27.19 

74 

32.63 

27.10 

74 

45.95 

27.22 

74 

33.14 


Table 4: Example 3; PSNR and CPU time for the best regularization parameter (maximum 


A = 0.04 
A = 0.1 
1/a = 0.0069 


FA-MD 

TV-MD 

FTVd 


.61 
29.87 
28.95 


15.95 

16.74 

0.73 


(f) Alg-Rectl (g) Alg-BC4 (h) Alg-BC4-NS (i) TV-MD (j) FTVd 


Figure 5: Example 4: true image, PSF, observed image, and restored images. 
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Algorithm 

PSNR 

Iter. 

CPU time(s) 

Regular, parameter 

Alg-BCl 

23.67 

51 

20.19 

a = 0.05 

Alg-BC3 

23.74 

17 

6.03 

a = 0.01 

Alg-Rectl 

23.59 

24 

12.44 

a = 0.02 

Alg-Rect2 

23.64 

17 

16.49 

a = 0.009 

Alg-BC4 

23.76 

17 

5.60 

a = 0.01 

Alg-BC4-NS 

23.53 

39 

14.78 

O 

II 

FA-MD 

23.44 


14.10 

A = 0.01 

TV-MD 

23.55 


13.37 

A = 0.11 

FTVd 

23.10 


0.89 

1/a = 0.0088 


Table 5: Example 4: PSNR and CPU time for the best regularization parameter (maximum 
PSNR). For our algorithms /x = 40. 



(a) Alg-BCl 


(b) Alg-BC4 


(c) Alg-BC4-NS 


(d) TV-MD 



(e) FTVd 





Figure 6: Example 4: 


Nort-East corner of the residual images. 


be computed by two antireflective transforms without requiring the PCG: in fact, the use of 
preconditioning would represent an unnecessary increase of the CPU time without increasing 
the PSNR with respect to Alg-BC3. Finally, Alg-BC4 is implemented by DCT like Alg-BCl. 

Table[5]shows that all the compared methods in this example provide comparable results. 
Nevertheless, it is interesting to observe that Alg-BC4 gives a slightly better restoration 
with a lower CPU time than the standard MLBA, i.e., Alg-BC3. TV-MD computes again a 
comparable restoration, but with more than a double CPU time. Figure [5] shows the restored 
images. 

To test the ability of the different algorithms in dealing with the boundary effects, FigurelB] 
shows the North-East corner of the residual images. We can see that Alg-BCl and FTVd 
have some ringing effects at the boundary, while Algorithm 4 and TV-MD do not show any 
particular distortion at the boundary. 


7 Conclusion 

In this paper, we have investigated several regularization preconditioning strategies for the 
MLBA applied to the synthesis approach with accurate restoration models for image deblur¬ 
ring and unknown boundaries. Our numerical results shows that Alg-BC4, which combines 
the favorite BCs (depending on the problem) with an approximated Tikhonov regulariza¬ 
tion preconditioner, represents a robust and effective algorithm. Indeed, it provides accurate 
restorations in all our examples with a reduced CPU time also in comparison to the state of 
the art algorithms nia and the standard MLBA when available, cf. Example 4. 
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We have investigated only the synthesis approach, but the same preconditioning strategies 
can be applied to the analysis and the balanced approach [61137] as well. Moreover, possible 
future investigations could consider the use of a preconditioner obtained by a small rank ap¬ 
proximation of the PSF as in |29] . strategies to estimate the parameter a, and nonstationary 
sequences to approximate the best a avoiding its estimation like in |28j . 
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